function merged = runfamamacbethreturns_quarterlyreturns_slope(merged,programFiles,Outfile,qtrdummy,allvarscalculated)


%prepare the variable for output file
outputFormat.Intercept = NaN*ones(500,1);
outputFormat.beta_lag_1q = NaN*ones(500,1);
outputFormat.fill_size_lag = NaN*ones(500,1);
outputFormat.fill_bm_lag = NaN*ones(500,1);
outputFormat.fill_lev_lag = NaN*ones(500,1);
outputFormat.fill_hhi_lag = NaN*ones(500,1);

outputFormat.fill_ol_lag = NaN*ones(500,1);
outputFormat.hhi_ol = NaN*ones(500,1);

outputFormat.fill_op_pft_1_lag = NaN*ones(500,1);
outputFormat.hhi_op_pft_1_lag = NaN*ones(500,1);

outputFormat.fill_hhi_slope_lag = NaN*ones(500,1);

%prepare interactions
merged.hhi_ol = merged.fill_hhi_lag.*merged.fill_ol_lag;
merged.hhi_op_pft_1_lag = merged.fill_hhi_lag.*merged.fill_op_pft_1_lag;
merged.fill_hhi_slope_lag = merged.fill_hhi_lag.*merged.high_slope;


outputFormat = struct2table(outputFormat);


%returns are already in percentages i.e., 8% is 8 in the file,
if qtrdummy == 1
    DepVar = 'qtrret'
    disp('Note: Running the regression every quarter')
    loc = merged.qtrdummy == 1;

    %     loc = merged.month_qtr == 3 | merged.month_qtr == 6 | merged.month_qtr == 9 | merged.month_qtr == 12;
    merged = merged(loc,:);

end



%rearrange data
merged = sortrows(merged,{'GVKEY','datadate'});


%prepare all specifications
Specifications{1} = {'beta_lag_1q','fill_size_lag','fill_bm_lag','fill_lev_lag'};
Specifications{2} = {'beta_lag_1q','fill_size_lag','fill_bm_lag','fill_lev_lag','fill_hhi_lag'};
Specifications{3} = {'beta_lag_1q','fill_size_lag','fill_bm_lag','fill_lev_lag','fill_hhi_lag','fill_ol_lag','hhi_ol'};
Specifications{4} = {'beta_lag_1q','fill_size_lag','fill_bm_lag','fill_lev_lag','fill_hhi_lag','fill_op_pft_1_lag','hhi_op_pft_1_lag'};
Specifications{5} = {'beta_lag_1q','fill_size_lag','fill_bm_lag','fill_lev_lag','fill_hhi_lag','fill_hhi_slope_lag'};

forregression(Specifications,DepVar,[],merged,outputFormat,qtrdummy,programFiles,Outfile);





function forregression(Specifications,DepVar,nlag,merged,results,qtrdummy,programFiles,Outfile)
if qtrdummy == 1
    nlag = 2;
    %for quarterly data, run the fama-macbeth using only quarterly

    disp('Note: Running the regression every quarter')
    disp(strcat('Lags: ', num2str(nlag)))
    loc = merged.qtrdummy == 1;

    merged = merged(loc,:);

end
counter = 1;
resultsMean = results;
resultsSD = results;
for j = 1:length(Specifications)
    disp(j)


    [out, nfirms, Mn, forMnpool, sdind] = performregressions(table2array(merged(:,DepVar)),merged,Specifications{j},nlag,merged.datadate);
    [~,~,ir] = intersect({'Intercept'},results.Properties.VariableNames);
    St = counter;
    En = counter + 2;
    results(St:En,ir) = array2table(out(:,1));
    sdind = sdind(:,3:end);

    for k = 1:length(Specifications{j})
        [~,~,ir] = intersect(Specifications{j}{k},results.Properties.VariableNames);

        results(St:En,ir) = array2table(out(:,k+1));
        resultsMean(St:St,ir) = array2table(mean(Mn(:,k)));
        resultsSD(St:St,ir) = array2table(mean(sdind(:,k)));
        resultsMean((St+1):(St+1),ir) = array2table(mean(forMnpool(:,k)));

    end
    disp(nanmean(nfirms))
    counter = En + 1;

    clear Mn forMnpool

end

%write the results in excel file
check = exist(strcat(programFiles,'\output\'),'dir');
if check == 0
    mkdir(strcat(programFiles,'\output\'))
end
cd(strcat(programFiles,'\output\'))
writetable(results,Outfile,'Sheet',DepVar);
writetable(resultsMean,Outfile,'Sheet',DepVar,'Range',strcat('A',num2str(3*(length(Specifications)+1))));
writetable(resultsSD,Outfile,'Sheet',DepVar,'Range',strcat('A',num2str(2*3*(length(Specifications)+1))));
cd(programFiles)

function [Final, Nfirms, Mn, forMnpool, sdind] = performregressions(y,xTable,Spec,nlag,dates)
%Inputs
%y: the dependent variable
%xTable: contains all the independent variables
%Spec: the specification for which we want to run the regression
%nlag: number of lags for newey-west

x = [];
for j = 1:length(Spec)
    disp(Spec{j})
    [~,~,ic]= intersect(Spec{j},xTable.Properties.VariableNames);
    disp(xTable.Properties.VariableNames{ic})
    x = [x, table2array(xTable(:,ic))];
end

x = [dates,y,ones(size(y)),x];
[ir,~] = find(isnan(x) == 1);
x(ir,:) = [];

dates = unique(dates);
forMnpool = [];
for j = 1:length(dates)

    loc = find(x(:,1) == dates(j));
    yy = x(loc,2);
    xx = x(loc,3:end);

    tx = [yy,xx];
    [ir,~] = find(isnan(tx) == 1);
    tx(ir,:) = [];
    [ir,~]  = find(isinf(tx) == 1);
    tx(ir,:) =[];


    yy = tx(:,1);
    xx = tx(:,2:end);

    R = rank(xx);

    Mn(j,:) = mean(tx(:,3:end));
    forMnpool = [forMnpool; tx(:,3:end)];


    if R == size(xx,2) && size(tx,1) > 7
        %         disp(size(yy))
        %         disp(size(xx))

        [b,~,~,~,tR2] = regress(yy,xx);
        R2(j,:) = tR2(1);
        Coeff(j,:) = [dates(j), b'];
        sdind(j,:) = [dates(j), nanstd(xx)];
    else
        if j == 1
            Coeff(j,:) = [dates(j), ones(1,size(xx,2))*NaN];
            sdind(j,:) = [dates(j), ones(1,size(xx,2))*NaN];
        else
            Coeff(j,:) = [dates(j) Coeff(j-1,2:end)*NaN];
            sdind(j,:) = [dates(j), sdind(j-1,2:end)*NaN];
        end
    end

    Nfirms(j,:) = size(tx,1);


end

disp(size(Coeff))
disp(sum(isnan(Coeff)))
disp(min(Coeff(:,1)))
disp(max(Coeff(:,1)))

Dates = Coeff(:,1);
Coeff = Coeff(:,2:end);
CoeffAll = nanmean(Coeff);
disp('number of firms')
nanmean(Nfirms)



%compute the newey west t-statistics
for j = 1:size(Coeff,2)
    y1 = [ones(size(Coeff(:,j))), Coeff(:,j)];
    [ir,~] = find(isnan(y1) == 1);
    y1(ir,:) = [];
    nwse = neweywest(y1(:,2) - nanmean(y1(:,2)),[],nlag);
    temptStat(1,j) = nanmean(y1(:,2))./nwse;
    tempSE(1,j) = nwse;
end

tStatAll_1 = temptStat; clear temptStat;
Coeff = CoeffAll; Signif = tStatAll_1; SE = tempSE;


% clear CoeffAll tStatAll;
Final = [];
for j = 1:size(Coeff,1)

    %     temp = [Coeff(j,:); Signif(j,:); Signif(j,:)*NaN];
    temp = [Coeff(j,:); Signif(j,:); SE(j,:)];
    Final = [Final; temp];

end



